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FCMtEWDRD 


This report presents results of the expansion and Improvement of the 
FORMA system for response and load analysis. The .acronym FORMA stands 
for FORTRAN Matrix Analysis. The study, performed from 16 May 1975 
through 17 May 1976 was conducted by the Analytical Mechanics Departoient, 
Martin Marietta Corporation, Denver Division, under the contract MASS* 
31376. The program was administered by the National Aeronautics and 
Space Administration, George C. Marshall Space Flight Center, Huntsville, 
Alabama under the direction of Dr. John R. Admire, Structural Dynamics 
Division, Systems Dynamics Laboratory. 

This report is published in seven volumes: 

Volume I - Programming Manual, 

Volume llA - Listings, Dense FORMA Subroutines, 

Volume IIB - Listings, Sparse FORMA Subroutines, 

Volume lie - Listings, Finite Element FORMA Subroutines, 

Volume lllA - Explanations, Dense FORMA Subroutines, 

Volume lllB - Explanations, Sparse FORMA Subroutines, and 
Volume me - Explanations, Finite Element FORMA Subroutines. 
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ABSTRACT 


This report presents techniques for the solution of structural 
dynamic systems on an electronic digital computer using FORMA (FORT RAN 
^trix ^alysis) . 

FORMA is a library of subroutines coded in FORTRAN IV for the effi- 
cient solution of structural dynamics problems. These subroutines are 
in the form of building blocks that can be put together to solve a large 
variety of structural dynamics problems. The obvious advantage of the 
building block approach is that programming and checkout time are limi- 
ted to that required for putting the blocks together in the proper order. 

The FORMA method has advantageous features such as: 

1. subroutines in the library have been used extensively for many 
years and as a result are well checked out and debugged; 

2. method will work on any computer with a FORTRAN IV compiler; 

3. incorporation of new subroutines is no problem; 

4. basic FORTRAN statements may be used to give extreme flexi- 
bility in writing a program. 

Two programming techniques are used in FORMA: dense and sparse. 
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LIST OF SamOLS 


[ ] 
{ } 

{ r 


T 


matrix 

column matrix . 

2 vector 
row matrix ^ 

transpose (%fhen symbol is a superscript) 


C ] 




®ij 


m designates the row size of matrix 
n designates the column size of suitrix 

a designates an element of oiatrix [aJ 
i designates the ith row of matrix [a] 

J designates the 1 th column of matrix [a] 



I. INTRODUCTION 


This volume presents an explanation of the function of each finite 
element subroutine in the FORMA library. Example problems are given in 
some cases to clarify the operations performed by a subroutine. 
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II. SOBBOirniiE EPLAimiOMS 


The subroutines are given in alphabetical order with nuwbers 
coming before letters. 



AXIAL 

Subroutine AXIAL calculates (on option) finite element: (1) mass 

matrices; (2) stiffness matrices (same as global load transformation 
matrices); (3) local load transformation matrices; (4) stress 
transformation matrices; and (5) vectors to locate the DOF (de~ 
grees of freedom) of the above matrices in the global DOF; for 
axial rod elements. The above matrices and vectors are written 
on disk units and constitute the output from this subroutine. 

All matrices are In dense programming logic. 

Each mass and stiffness matrix, size 6x6, is in the global co- 
orlnate directions. The global coordinate order for each element 
is (U,V,W) joint 1, then joint 2 where U, V, W are translations. 

If the Euler angles are zero at a joint, then U * 6„, V = 6 , 

w - X Y 

Each global load transformation matrix, size 6x6, relates loads 
at the jod ends In the global coordinate directions to deflections 
in the global coordinate directions. The row order in this 
matrix is (Py, P^, P^) joint 1, then joint 2 where P is force. 

Each local load transformatipn matrix, size 2x6, relates loads at 
the rod ends in the local coordinate system to deflections in the 
global coordinate directions. The row order* in this matrix is 
^xi’ ^x2 axial force. 

Each stress stransformatlon matrix, : ze 2x6, relates stresses at 

the rod ends In the local coordinate system to deflections in 

the global coordinate directions. The row order in this matrix 

Is 0 , , a „ where a is normal stress, 
xl x2 

Each location vector (IVEC) locates the DOF of each finite element 
in the global DOF. For example, IVEC(6)»834 places element DOF 6 
into global DOF 834. IVEC(3)=0 omits element DOF 3 from global 

DOF. This constrains element DOF 3 to zero motion. 

The above matrices are calculated by using joint data and element 
data. The joint data is obtained from three matrices input to 
this subroutine which are (1) joint global X, Y, Z locations; 

(2) joint global DOF numbers; and (3) joint Euler angles. 

The element data read in this subroutine is (1) options for mass, 
stiffness, local load transformations, stress transformations; 

(2) element material properties; and (3) element joint numbers, 
cross-sectional area. 

Each mass matrix is calculated by transfer to subroutine MASIA. 

Each stiffness matrix, loads, and stress transformation matrix 
is calculated by transfer to subroutine STFIA. 



lAl 


Subroutine BlAl calculf tes a buckling (sonetlmes referred to as 
geometrical stiffness, initial stress, or stability) matrix for 
an axial rod element with unrestrained boundaries. The buckling 
matrix is based on a unit axial load. The buckling matrix is in 
the local coordinate system of the rod. 

DE SCRIPTION OF TECHNIQUE 

From Theopy cf r!atp-ix StPuatuml Analysis by J. S. Przemienieckl, 
McGraw-Hill 1968, we obtained the buckling matrix. The strain 
energy for buckling is obtained as 



where the kernel matrix is the buckling matrix. A unit axial load 
of F = 1 is assumed here. L is the rod length. 

The degrees of freedom are shown in the following sketch. 


z 




BU2 


Subroutine B1A2 calculates a buckling (sometimes referred to as 
geometrical stiffness, initial stress, or stability) matrix for 
a beam element %rith unrestrained boundaries. The buckling matrix 
is based on a unit axial load. The buckling matrix is in the 
local coordinate system of the beam. 


DESCRIPTION OF TECHNIQUE 


From Theoi^y of Matrix Struaticral Analysis by J. S. Przemleniecki, 
McGraw-Hill 1968, we obtained the buckling matrix. The strain 
energy for buckling is obtained as 




where the kernel matrix is the buckling matrix. A unit axial 
load of F«1 is assumed here. L is the beam length. 


The degrees of freedom are shown in the following sketch. 




BAR - 1/2 


Subroutine BAR calculates (on option) finite element (1) mass 
matrix; (2) stiffness matrices (same as global load transformation 
matrices); (3) buckling matrices for unit load; (4) local load 
transforma tlort matrices; (5) stress transformation matrices; and 
(6) vectors to locate the DOF (d#^ ,rees of freedom) of the above 
matrices in the global DOF, for . mbined axial-torslon-bending 
bar elements. The above matrices and vectors are written on disk 
units and constitute the output from this subroutine. All matrices 
are in dense programming logic. 

Each mass, stiffness, an » buckling matrix, size 12x12, is in the 
global coordinate directions. The global coordinate order for 
each element is (U, V, W, P, Q, R) joint 1, then joint 2 where 
U, V, are translations and P, Q, R are rotations. If the 
Euler angles are zero at a joint, then U * 6 , V = 6 , W » 6 , 
p = Q = 6^. » = 6^. 

Each global load transformation matrix, size 12x12, relates loads 
at the bar ends in the global coordinate directions to deflections 
in the global coordinate directions. The row order in this matrix 
is (Py, P^, P^, Mp, M^, Mj^) joint 1, then joint 2 where P is force 

and M is moment. 


Each local load tran£ lormation matrix, size 12x12, relates loads at 
the bar ends in the local coordinate system to deflections in the 
global coordinate directions. The row order in this matrix is 


P P 
xl* x2* 


M 


xl’ ^2’ ^yl’ ‘y2’ 
P is force and M is moment. 


M 


M 


M 


zl* ' z2* zl’ z2* ‘yl* 


M where 
y2’ 


Each stress transformation matrix, size 12x12, relates stresses 
:'.t the bar ends in the local coordinate system to deflections in the 
global coordinate directions. The row order in this matrix is 



^x2 

M * r. 

xl 1 

\2* 

h 

A, ' 

A. * 

J 



1 
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1 

2 



!z2 

M * c , 

M^2* 


A, ’ 

A^ • 

1 , 

1 ^ 


1 

2 

zl 

z2 



^.2 

M * c , 
yl zl 


c ^ 
zz 


^2 ’ 
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BAR ~ 2/2 

vkmtm F is fores, M is wissot, A is cross-soetioasl. sros, r is 
Bistsscs froa torsi<m sxis ^ ooCsr fiber, J Is cross* 
ssctlon SaiBC VMSt's Corsioa consrsnt is c is disumes 
froa bending nsstrsl gins to ostsr fibsr, sod ^ is srss 
■oa^ df iasreis. 

B^i locscion vsetor (IFEC) locstss die VOt of fislts sl s as o t 
is Che global BCff. For ssaopls, OTC<b)«834 glsess slsasoc DOF b 
into global BOF B34. IVBCC3)«0 sales slsaiat S3F 3 froa global 
IXMP. This oooscrains e l sas s t DW 3 .to ssro astisa. 

These aacricss are caloilaCed by using Joint data and elcaent 
data. The joint data is obtained froa three aatrices input to this 
storoutine whl<h are (1) Joiac global X, T, Z locacions; (2) joint 
global DOF nadisrs; and (3) joint Euler angles. 

The eleasnt data are (1) tobi<^ for aass, stiffaess, local load 
cransforaatiems, stress transforaatioos; (2) eleaent aaterlal 
properties, and (3) eleaent Joint mabers. 

Each ness aatrix is calculated by transfer to subroutine HASH. 

Each stiffn^s aatrix, loaito, and stress transforaation aatrix is 
calculated by transfer to subnmtine STFIB. Each unit load buckling 
aatrix is calculated by transfer to subroutine BOCIB. 



BOCIB > 1/3 

Subroutia* BOCIB calculatu • finite eleamit buckling (sonetiMS 
referred to as geoaetrical stiffness, initial stress, or stability) 
■atrix for a coabbied axial-torsion<-bendlng bar elMmt with un- 
restrained boui^aries. The buckling matrix is based on unit axial 
' load. 

The buckling aatrix, sise 12x12, is in the global coordinate directions. 
The global coordinate order for each element is (D,V,W,P,Q,g) Joint 1, 
then joint 2 where U,V,U are translations and P,Q,R are rotations. If 
thm Euler angles are sero et a joint, the 

Q"0y. 

This aatrix is calculated by first computing a buckling aatrix in 
the local coordinate system for either an axial rod (where the 
buckling matrix Is sometimes referred to as the string stiffness) 
or a beam. A direction cosine aatrix is then used to transform 
the buckling matrix from the local coordinate system to the global 
coordinate system. 

DESCRlPTKMi OF TECHHKXJE 

The calculation of the buckling aatrix in the global coordinate 
directions is accomplished as follows. First a buckling matrix 
is calculated in the local coordinate system for either an axial 
rod (reference Subroutine BlAl) or a beam (reference Subroutine 
B1A2). 

k sketch of the bar is given for reference as 





X 



BUCIB - 2/3 

Til* atraln ancrgy for buckling or goomtric •tiffiom* using local 
coordinstM is 

D - I [bj^l i\} UJ 


uhsro 



b^^ refers to terms from Subroutines BlAl or B1A2. 

The deflections in the local system are related to deflections In the 
global systoD by 


{hj^} • [X] {hg} 


12 ] 



where [y] Is e direction cosine ■etrlz (reference Subroutine 
DCOSIB) Including Euler angles, slse 12x12, and 

(h„)^ - \ «i »2 “2 '2 "2 *2>- 

U, V, W are translations and P, Q, R are rotations. 
Substituting Eq [2] Into Eq [1] gives 

“ ■ 1 '"ol "«c> 

T 

where [b^] “ [y] [b^^] [y] Is the buckling matrix In' global 

' coordinate directions . 



DCOSIA > 1/2 

Subroutine DCOSIA calculates a direction cosine natrix for an 
axial rod elenent. This matrix relates local coordinate displace- 
ments to global coordinate dlsplacoaents. Euler angles at each 
joint are included. Global X, Y, Z coordinates and Kiler angles 
at each of the two rod ends are needed for this calculation. 

DESCRIPTION OF TECHNI(»E 

A sketch of the rod is given for reference as 


X (origin at point 1) 



The vector P between points 1 and 2 is 

f r . p, j . Pj K 

The unit vector is then 

- [Pj^ T + P^ J + P^ K]/)lp 
where 



11 } 

[ 2 ] 





DCOSU - 2/2 


The coefficients of J, K for the unit vector Cp ere the 
direction cosines of the line between points I and 2. 

The relation between local and gldbal X, Y, Z displaceamts is 
then 


xl 
where 


tep] 


0 


^^XY2^2 


13) 




A 3x3 Euler angle transformation matrix (reference subroutine EULER) 
relates global X, Y, Z translations to global U, V, W translations 
at each joint. That is. 


— - 
^"^XYZ^l 


■[Ell 

0 


{UVWlj 

kPb M 


0 

[Elj 


1 

a 


where 


{UVW} 


Substituting Eq [4) into Eq [3] gives the direction cosine matrix. 


(A) 


1 

o> 

X 

1 


0 ■ 


"{UVW}j^' 


- 



{UVW >2 


[ 5 ] 


where 



DCOSIB - 1/5 


Subroutine DCOSIB calculates a direction cosine matrix for a 
combined axial-torsion-bending bar element* This matrix relates 
local coordinate displacements to global coordinate displacements. 
Euler angles at each joint are included. Global Y, Z coordinates 
and Euler angles at each of the two bar ends plus coordinates 
of a reference point are needed for this calculation. The ref« 
erence point defines the local xy plane. 

DESCRIPTION OF TECHNIQUE 

A sketch of the bar is given for reference as 



The vector P between points 1 and 2 is 

P - K (1) 

The unit vector is then 


[p^ t + p^. j + p. tj/.p 


121 



DCOSIB - 2/5 


where 


P 


X 


X 


1 



^2-^1 


'p -V^X ♦ "v + ‘■z 

The coefficients of 1, T, K for the unit vector*^ are thv- direction 
cosines of the line between points I and 2. 

The vector cross poduct P x 1, 3 will give a vector (R)J_ plane 
1. 2, 3. 

R = ^ 1^3 = iPjj T + Py J" + P 2 Kl l(X 3 -X^) 1 + (Y 3 -Y 3 ) J + ( 23 - 2 ^) t) 

-B^T^R^r+R^K. 

The unit vector along R is then 

= [R^T+ Ry 7+ R^ll/lj^ [4] 


where 


“* ■ h «3-^> - ’’Z 
^ «3-*l> - ‘’x <^3'^1> 

’>Z • ‘'x ''^3-'^l> - ’'v (*3-*l' 



The coefficients of I, J, K for the unit vector are the direction 
cosines of a lir^J^ plane 1, 2, 3. 

The vector cross product R x P will give a vector (Q) JL line 1, 2 in 
the plane 1, 2, 3* 


Q = R > P » IR^ T 4 R^ J + R^ K] X 

Q^Kl . 


[51 



The unit vector along Q is then 
where 


DCOSIB - 3/5 


[ 6 ] 


" ®X ^Z 

Qz - "y - \ "X 

t, -VqI + q§ + 

The coefficients of 1« J, K for the unit vector e^ are the direction 
cosines of a_l_line 1, 2 and in the plane 1, 2, 3. 

The relation between local and global X» Y, Z displacements is then 




■- — 
[ep] 


^"^XYZ^l 

CM 

<o 


[epl 


^®XYZ^1 

«xl 


tep] 



®x2 


lap] 


^®XYZ^2 

^1 


tep] 


^2 


lepl 


®zl 




®z2 




/zl 




'z2 


1«rJ 




[eg] 


1 

a> 

m: 

K) 


t^Qi 




where 



and {6 xyz> " 

1 

CD 

X 

6 



Y 


Y 

6 



Z 


Z 


DCOSIB - 4/5 

f7l 


A 3x3 Euler angle transformation matrix (reference subroutine EULER) 
relates global X,Y,Z displacements to global U,V,W translations and 
P,Q,R rotations. That Is, 


p'^XYZ^l 


■(Ell 


{UVWlj^ 

CD 


(Ell 


{PQRlj^ 

^*XYZ^2 

K 

(EI 2 


{UVW}^ 

_^®XYZ^2_ 


(El2 


{PQR}2 


where {UVW} * 

U 

and {PQR> = 



V 


Q 


W 


R 


Substituting Eq [8] Into Eq [7] gives the direction cosine matrix. 


» — 
^<1 


■(e,li 


"(UVWlj^" 

^2 


(«pl2 


{PQRlj^ 

®xl 


[epli 


{UVW}^ 

^2 


(epl2 


{PQR}2 



f^qll 


^2 


t«Ql2 

♦ 


®zl 


f«Rll 


^z2 




^1 


f^RJl 


'z2 


^®R^2 


°yi 




1 

04 

Ds 

0 

1 


teQl2_ 



[ 8 ] 



where 

I«plj • {«pl CElj «tc. 


DCOSU - S/5 



DC0S2 > 1/6 

Subroutln* DC0S2 calculatas a diracclon cosina natrlx for a caa~ 
bined membrane-bending triangle plate element. This matrix relates 
local coordinate displacements to global coordinate displacements. 

Euler angles at each joint are Included. Globfd X, Y. Z coordinates 
and Euler angles at each of the three comers are needed for this 
calculation. 

DESCRIPTlCai OF TECHNIQUE 

A sketch of the triangle plate is given for reference as 


3 



The vector P between points 1 and 2 is 

P - T + Py J + K tl] 

The unit vector is then 


[ 2 ) 



DC0S2 - 2/6 


where 


* *2 ” 



The coefficients of I, J, K for the unit vector Cp are the 
direction cosines of the line between points 1 and 2. 

The vector cross product P ^ 1, "i will give a vector (R)J. plane 

1, a. 3. 

R - 7 X ITT - [Px*! + ^ ^ 

. + Rj T + Rj K. ,3, 

The unit vector along R is then 

tjj - iRjjt + Ry"! + Rj^T]/£p [4] 

where 

“x - 'y \ 

»Y ■ fz <V*1> - 'X 

“z ■ ■’x 



The coefficients of I, J, K for the unit vector are the direction 
cosines of a llnej_plane 1, 2, 3. 



DC0S2 - 3/6 

The vector cross product R P will give a vector W)^linv 1, 2 
In the plane 1, 2, 3. 


Q-R»P-lRjjI + ■■ *2^1 " 'y ^2 

■Q^t+q^T+Q^T. 

The unit vector along R is then 
where 

Ox ■ 'y "z - ■Si ^ 

^Z 

^Z ” 

‘q -l/Ox + Oy^ + Q| . 


15] 


16] 


The coefficients of I, J, K for the unit vector e^ are the direction 
cosines of a line line 1,2 and in the plane 1, 2, 3. 



DC0S2 - 4/* 

The r«latlon b« tw « n local and gl^»al Z» Y» Z dlaplaoaanta Is 
then 






V 






Y 


Y 




^ 2 ^ 


^ 2 ^ 


DC0S2 - 5/6 


A 3x3 Eul«r angle transformation matrix (reference subroutine EULER) 
relates global X, Y, Z displacomnts to global U, V. W translations 
and P, Q, R rotations. That is. 


’t*XY 2 >l 


"(Elj 


{OVW>^ 

t«m>i 


[Ell 


{Pcpilj 

^^XYZ ^2 

m 

[Elj 


(UVWlj 

^®XYZ ^2 


[Elj 


{PQR >2 

^*XYZ^3 


[EI 3 


{UVWlj 

(e } 

' XYZ^3 


l=l3_ 


{PQRlj 


where (UVU) - 

• « 
U 

and 

{PQR} - 



V 



Q 


w 

• m 





[ 8 ] 


Substituting Eq [ 8 ] into Eq [7] gives the direction cosine matrix. 














EDL^ - l/A 

Subroutine EULER celculaces the Euler rotation trensforiiation 
Dtatrlx such that 


*] 


u 


- [T1 j 

V 

. Z ] 


w 


where 


X 
Y 
. Z 


global X, Y, Z coordinate system. 


' U 
V 

w 


rotated U, V, V coordinate system. 


[T] B Euler rotation transformation based on a global d , 6 , 
and permutation. 

DESCRIPTION OF TECHNIQUE 

The first Euler rotation is 0 about X to form the X', Y', Z' 
coordinate system. 



Y 



BOLOL - 2/4 


The relaticm^lp betweai the two coordinate systeas can he 
written as 



IT6,1 



%mere 


fi 0 


0 




0 

0 cos0^ 


The second Euler rotation is 6^ about to fom the Z 

coordinate systoi* 



The relationship between the two coordinate system can be written 
as 




EDLE& - 3/4 


where 


cosS. 


sine.. 


tte^l 


-sine. 


1 

0 


co*e„ 


The third Euler rotetion is 6 about Z" to fora the U, V, U 
coordinate systen. 



The relationship between Che two coordinate systens can be written 
as 


(*"> 1 

1“ 

Y"}- IHjI 


( 2"" ) ^ 

Iw) 


where 

|T6jl 


cos6 -sin6 0 

«# i» 

sinS^ cose^ 0 


0 


0 


1 



EDLER - 4/4 

The conplat* Euler rotation transforaetion can be written as 
IT) - ixejj] iTe^l iTejjl 

Per€ondng the three multiplicatlona reaulta in 



cosd^ 

coae^ 


-0080^ 

8in0z 


sinO^ 



aine^ 


“•*x 

00802 



OO80y 

IT) • 

+sin®jj 

ainS^ 

00802 


sinO^ 

aln02 





aine^ 


*‘“'x 

00802 


co.6jj 

c..«. 


-“•®x 

aine^ 

00802 

+00802 

sin02 

8in0z 





FINEL 


Subroutine FINEL calculates (on option) using finite elements: 

1) an assembled mass matrix; 2) an assembled stiffness matrix; 

3) element local load transformation matrices; 4) element glo* 
bal load transformation matrices; !>) element stress transforma- 
tion matrices; 6) element unit load buckling matrices; and 7) 
vectors (IVEC) to locate the DOF (degree of freedexa) of the 
element matrices in the global DOF. 

The types of finite element available (and the related 
subroutine) are axial rod (AXLAL) , combined axial- torsion- 
bending bar (BAR)^ triangular plate (TRNGL), quadrilateral 
plate (QUAD), rectangular shear panel (RECTSP), tetrahedron 
(TETRA), and pentahedron (PENTA) . Tlie subroutine to be used 
is specified by reading this information (e.g., AXIAL, BAR, 
etc.) from an input data card. 

The assembled mass and stiffness matrices are output from 
this subroutine in sparse FORMA subroutine format on disk units. 
The DOF order is specified by a joint degree- of- freedom matrix, 
[\JD0F] , which is input to this subroutine. 

Tlie element matrices and vectors are in dense programming 
logic and written on disk units as output fromthis subroutine 
also. The sizes of these element matrices and vectors are de* 
termined by the specific finite element used. Each vector 
(IVEC) locates the DOF of each finite element in the global 
DOF. For example, IVEC(6)*834 places element DOF 6 into global 
DOF 834. IV^C(3)*0 omits element DOF 3 from global DOF. This 

constrains element DOF 3 to zero motion. 

The finite element matrices are calculated by using joint 
data and element data. The joint data, obtained from three 
matrices input to this subroutine, are 1) joint global X, Y, Z 
locations; 2) joint global DOF numbers; and 3) joint Euler 
angles. The element data is read in the specified finite ele- 
ment subroutine. Reference AXIAL, BAR, etc. for this data. 

Assembly of the element mass (or stiffness) matrices into 
the assembled mass (or stiffness) matiix for the total structure 
is accomplished by FORMA subroutine YRVAD2 to obtain the sparse 
subroutine format. 



KUl - 1/4 

Subroutine KlAl calculates a stiffness matrix and stress trims-* 
formation matrix for an axial rod element with unrestrained 
boundaries* The stiffness matrix is in the local coordinate sys- 
tem of the rod. The elements of the stiffness matrix represent 
the distributed stiffness properties of the rod. These elements 
are calculated by assuming constant axial force. The stress 
transformation matrix relates stress at the rod ends in the local 
coordinate system to deflections in the local coordinate system. 

The rod may be linearly tapered or uniform. 

DESCRIPTION OF TECHNIQUE 

The replacement of the distributed axial stiffness of a rod by 
a stiffness matrix is obtained using a strain energy approach as 
follows . 

Consider a rod that is loaded with an axial force at point 1 
and restrained at point 2 as shown in the sketch. 



X 


The strain energy is defined by 


U 



P(x)^ 

A(x) eU) 


dx 


Where 


fii 


P is the axial force, 

A is the cross-sectional area, 

E is Young's modulus of exasticity 

X is the local coordinate system and longitudinal axis of the rod. 
The origin is at point 1, that is, x^^ * 0; X 2 ® L (rod length). 
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To Integrate Eq [1], the axial force is assuaed constant and equal 
to the axial force at point 1, that is, 

P(X) - ?y (21 

Young's modulus of elasticity is also assumed constant, that is, 

E(x) - E. (31 

The cross-sectional area is assumed to vary linearly, that is, 

A(x) * + X (A^-Aj^)/!. (41 


Substituting Eq [21 through [41 into Eq [1] gives the strain 
energy as 


U 


1 lit 

2 E 1 


Aj^ + X (A2-A^)/L 


2 (A^-Aj^) E 


dx 


(51 

(61 


Application of Castigliano's theorem gives the axial deflection 
of point 1 relative to point 2 as 


A6 = 


3U 

3P, 




from which 

(A2-A^) E 

^1 “ L tn (Aj/A^) 


A6 


(71 


The restraint at point 2 is removed by application of the trans- 
formation 


A6 - [1 - 1 ) 


(81 
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where 5^ and displacements along the rod x-axls at 

rod*ends 1 and 2, respectively. Substitution of Eq [7] and 
Eq [8] Into Eq [6] gives the strain energy for a rod with un- 
restrained bounterles as 


U 


1 

2 


* 2 > 


*1.1 *1.2 


■‘x‘ 

(sym) Zj 2 




where 

(A^-Aj) E 

*1,1 “ L In (A 2 /A^) 


* 1,2 * ■ * 1.1 
* 2,2 “ * 1,1 


[9] 


1101 

111 ] 

tl2] 


The kernel matrix of Eq [9] is the stiffness matrix that re- 
presents the axial stiffness of a rod with unrestrained boundaries. 


For constant cross-sectional area, l.e. , A^ > A 2 


A, Eq [10] Is of 


indefinite form. For this case integration of Eq [5] yields 
1 


[6a 1 


from which 
AE 

*1,1 “ L * 
as before 


[10a] 


z 


1,2 


z 


1,1 


* 2,2 ” * 1,1 


The elements in the stress transformation matrix are easily calculated. 
The following sketch shows the sign convention. 



1 


2 
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The rod end forces (P^, F 2 ) can be expressed in terms of the 
rod end displacements (5^, as 



* 1.1 * 1,2 


"sym) z 


2.1 


This is obtained from Eq [7] or applying Castigliano's theorem to 
Eq [9]. The stress at the rod ends is simply 


[13] 


s 


1 


Pj^/A 


1 


and 


'2 ■ . 


or 


r n 

*1 

“ [T 1 



8 


^2 


m 


where 


tTj 


* 1 , 1^^1 * 1 , 2^^1 
* 2 , 1^^2 * 2 , 2^^2 


Is the stress transformation matrix. 


[14] 


[151 


8 ^. 8 ^ will be opposite in sign. Tension and compression in the 
rod Is determined as follows: 

Tension: Sj^ (-) , 82 (+) 


Compression: s^^ (+), 82 (-) 
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Subroutine KlBl calculates a stiffness matrix and stress trans- 
formation matrix for a bending (plus shear) beam element with 
unrestrained boundaries. The stiffness matrix is in the local co- 
ordinate system of the beam. The elements of the stiffness matrix 
represent the distributed stiffness properties of the beam. These 
elements are calculated by assuming uniform shear and linear bend- 
ing moment variation. The stress transformation matrix relates 
stress at the beam ends In the local coordinate system to deflec- 
tions in the local coordinate system. The bean may be tapered or 
uniform. 

DESCRIPTION OF TECHNIQUE 

The replacement of the distributed bending and shear stiffness of 
a b»am by a stiffness matiix Is obtained using a strain energy 
approach as follows: 

Consider a beam that Is loaded with a shear and moment at point I 
and restrained at point 2 as shown in the sketch 


z 




/ M(x)^ 

V(x)l \ 

1 . 

^E(x) I(x) 

KA(x) G(x) 1 


dx 


til 
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where 

H(x) is the bending moment, 

V(x) is the shear, 

E(x) is Young’s modulus of elasticity of the material, 

l(x) is the cross-sectional moment cf inertia about the 
beam's neutral axis, 

K is the shape factor (e.g., K ■ 1 for a solid circular 
cylinder, K * 0.5 for a thin walled circular cylinder), 

A(x) is the cross-sectional area, 

G(x) is the shear modulus of elasticity of the materiel, 
and 

X is the local coordinate system and undeformed lo igltuolnal 

axis of the beam. The origin is at point 1. That is x. ■ 0; 

X 2 ■ L (rod length) . 

To Integrate Eq [1], the following assumptions ere made. First, 

the shear is assumed oomtant and equal to the sV>eer force at point 

1, that is, 

V(x) - [2] 

Second, the bending moment is assumed to vai'y Zineax'Xy, that is, 

M(x) » X [3] 

Third . 

I(x) and A(x) are assumed to vary linearly, that is 


I(x) - + X (I 2 -Ij^)/L M 

and 

A(x) ■ Aj^ + X (A2~Aj^)/L t^b] 

Fourth, the moduli of elasticity, E and G, are assumed constant, 
that is 

E(x) - E [5a] 

G(x) - G [5b] 


Substituting Eq [2] through [5] Intv/ Eq [1] gives the strain energy 



0 - J iVj M^l j ^E{I^ + X 

l: :])* [::] 

• I l»l «1> '“1 pi] 

' 22 J Kl 


r:| 
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»C{A, X (A -A )/L} 


vhere 


11 El 


(R-1) 1^2 ■ R-1 ^ * KGCAj-A^) A^ 

^ (R-i) ~ ^ ^ 


12 El 


^22 “ EIj^ (R-1) ^ 

R - I^/l, 

For constant bending stiffness, i.e., El^ » EI 2 “ El, and constant 
shear stiffness, i.e., KA^G “ KA 2 G * KAG, Eq (7a] [7b] [7c] are 
of indefinite form. For this case, integration of Eq [6] yields 

fj^j^ * L^/3EI + L/KAG 

fj2 - L^/2EI 

f22 = L/EI 

Application of CastigHaao 's theorem to Eq [7] gives the lateral 
translation and rotation of point 1 relative to point 2 as 



[ 9 ] 
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Solving Eq [9] Cor and and substituting into Eq [7] gives 
the strain energy as 


C22/D 

(syn) [Ae 


1101 
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The kernel natrlx of Eq [12] is the stiffness natrix that repre- 
sents the bending and shear stiffness of a bean with unrestrained 
boundaries. 

The elements in the stress transformation matrix are calculated 
as follows. The following sketch shows the sign convention. 



Applying Castigliano’s theorem to Eq [12] gives the forces at the 
beam ends in terms of the displaconents at the beam ends, that is. 


Sli/ 3d 




X 

3U/36, 

St 


- [2] 


3U/39j^ 




®1 

3U/302 

_ _ 


«2 


®2 


The shear stress is calculated as 
T » V/A 

and the bending stress is calculated as 
0 - Mc/I 


where c is the distance from the beam's neutral axis to the outer 
fiber. 



Therefore, 




v^/*i 



^2 


CM 

> 


^2 


- 

MiCi/li 

- ITS] 

*3 

'’z 

-J 






where 


[TS] 


Z (row 1)/Aj^ 

Z (row Dtk^ 

Z (row 3) 

Z (row 4) c^/Ij 
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is the stress transformation matrix 
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Subroutine KlCl calculates a stiffness aatrlx and stress trsns- 
foraation matrix for a torsion rod elsnent with unrestrained 
boundaries. The stiffness matrix is in the local coordinate sys- 
tem of the rod. The eleioents of the stiffness matrix represent 
the distributed stiffness properties of the rod. These el^ients 
are calculated by assuaing constant torque. The stress trans- 
formation matrix relates stress at the rod ends in the local co- 
ordinate system to rotations in the local coordinate syst«. The 
rod may be tapered or uniform. 

DESCRIPTION OF TECHNIQUE 

The replaceaent of the distributed torsion stiffness of a rod by 
a stiffness matrix is obtained using a strain energy approach as 
follows. 

Consider a rod that is loaded with a torque T^ at point 1 and 
restrained at point 2 as shown in the sketch. 



Rod Element 


The strain energy is defined by 


U = 



T^(x) 
J(x) G(x) 


dx 


[ 1 ] 


where 

T is the torque 

J is Saint Venant’s torsion constant » 

J = 7rR**/2 for a solid circular section 
J * 2rrR^t for a thin walled circular section 

G is the shear modulus of elasticity 

X is the local coordinate system and longitudinal axis of the rod. 
The origin is a point 1, that is Xj^ * 0; X 2 ® L (rod length). 
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To integrate Eq (1], the torque is assueed constant and equal 
to the torque at point 1, that is, 

Kx) - [21 

The shear nodulus of elasticity is also assuned constant, that is, 

G(x) - G (31 

Saint Venant's torsion constant is assusted to vary linearly, that is. 


J(x) - + X 141 

Substituting Eq [21 through [4] into [1] gives the strain energy 
as 


V 


I r| j-t 

1 1 L 

2 ® 


1 

xUj-Jj^) /L 

in 


dx 


[51 

[61 


Application of Castigliano's theorem gives the rotation of point 1 
relative to point 2 as 


A6 





G 

from which T, ■ - — - — 7~r~7~r~\ I /I 

1 L in (J 2 /J^) 

The restraint at point 2 is removed by application of the 
transformation 


AS = [1 -1] 



[8] 


where 6^ and 6^ are the x~rotations at rod ends 1 and 2, respectively. 

Substitution of Eq [7] and [8] into [6] gives the strain energy for 
a rod with unrestrained boundaries as 


“ • I '»! *2l 


'l.l 'l,2 


(sym) z 


2,2j 


1 

’2 


[9] 



where 


(Jj-Jj) G 


*1.1 ■ 


2 » 2 

1.2 * 1,1 


' 2,2 " * 1.1 
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[101 

[111 

[121 


The kernel netrlx of Eq [9] is the stiffness natrix that repre- 
sents the torsional stiffness of a rod with unrestrained 
boundaries. 


For a constant cross section, i.e., * *^2 * 

indefinite form. For this case, integration of Eq [5] yields 

L 

U > — — 

2 JG 


(6al 


from which 



as before 


tlOa] 


z 


1,2 


z 


1,1 


* 2,2 " * 1,1 

The elements in the stress transformation matrix are easily calculated. 
The following sketch shows the sign convention. 


T 


1 ' 


1» 



X 


The rod end torques (T^^, T 2 ) can be expressed in term's of the rod 
end rotations (6^, as 


'^1 

IK 

*1,1 *1,2 


®1 

"2 

be m 


(sym) z^ 2 


«2 


[ 13 ] 
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This Is obUlnsd from Eq [7] or applying Cutlgllaho's theoram to 
Eq [9]. Tha maxi mum stress is In the outermost fiber and Is 


*1 ■ ■'i h'h 


and 

Sj - 


or 

fs 


t^el 


where 


II, 1 


‘i.i '■I'-’i 


* 2.1 * 2^'’2 


*1,2 V^l| 

*2,2 V'’2| 


Is the stress transformation matrix. 
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Subroutine K2A1 calculates a stiffness matrix and stress trans- 
formation matrix for a membrane triangle plate element with 
unrestrained boundaries. The stiffness matrix is in the local 
coordinate system of the triangle plate. The elements of the 
stiffness matrix represent the distributed stiffness properties 
of the triangle plate. These elements are calculated by assuming 
a quadratic displacement (linear strain) field. The stress trans- 
formation matrix relates stresses at the triangle vertices in the 
local coordinate system to deflections In the local coordinate 
system. 

DESCRIPTION OF TECHNIQUE 

The replacement of the distributed meodirane stiffness of a triangle 
plate by a stiffness matrix is described in DM 109 Linear Strain 
Membrane Triangle Element by W. A. Benfxeld and C. S. Bodley. 

The triangle is Illustrated L. the sketch with the degrees of 
freedom shown 


e 


zl 



X 


The order of the degrees of freedom is 


f[6 6 9 ], [6 <S 9 ], [6 6 6 ],1 . 

l x y z'l x y z 2 x y z 3\ 

The row order of the stress transformation matrix is 

[[n c X ] fo a r ). [o a x ]_1 

L x y xy 1 X y xy 2 x y xy 3j 


where a is normal stress and x is shear stress. 



K2B1 - 1/2 

Subroutine K2B1 calculates a stiffness matrix and stress trans- 
formation matrix for a bending triangle plate element with un- 
restrained boundaries. The stiffness matrix is In the local co- 
ordinate system of the triangle plate. The elements of the stiff- 
ness matrix represent the distributed stiffness properties of the 
triangle plate. These eleaents are calculated by assuming a cubic 
displacement (linear curvature) field. The stress transformation 
matrix relates stresses at the triangle vertices in the local 
coordinate system to deflections in the local coordinate system. 

DESCRIPTION OF TECHNIQUE 

The replacenent of the distributed bending stiffness of a triangle 
plate by a stiffness matrix uses the technique (essentially) 
described in Triangular Elements in Plate Bending by C. P. Bazely, 

Y. K. Cheung, B. M. Irons, and 0. C. Zlenkiewicz, AFFDL-TR-66-80 , 
November 1966. The triangle is illustrated in the sketch with the 
degrees of freedom shown 



The order of the degrees of freedom is 

f(6 9 0 ], 16 6 0 1, [6 0 e 1_] . 

z X y 1 z X y 2 z x y 3J 



K2B1 - 2/2 


The row order of the stress transforaetlon matrix is 


Ft° 0 T ], 

L X y xy 1 

Flo 0 T ], 
L X y xy 1 


I®x ^xy^2 ^'^x °y ^xy^sj end 

‘‘’x’yVz '°x''y Vs]" 


where 


0 is normal stress 
T is shear stress 
t is plate thickness. 



Subroutine MlAl calculates a lumped mass matrix for an axial rod 
elemant with unrestrained boundaries. The mass matrix is in the 
local coordinate system of the rod. The elements of the mass 
matrix represent the distributed mass properties of the rod. 

These matrix elements are calculated by lumping the rod's total 
mass to the rod end points using static equivalent forces. The 
rod may be linearly tapered or uniform. 

DESCRIPTION OF TECHNIQUE 

The replacement of the distributed mass of a rod by a mass matrix 
is obtained by lumping the rod's total mass to the rod end points. 
This lumping is equivalent to static beaming. The rod is illus- 
trated in the sketch in which the rod cross-sectional area varies 
linearly, that is, A(x) “ (A2~Aj^) . 



The total mass of the rod is 
M » p L (A^+A2>/2 

where 

p Is the mass density 

A Is the cross-sectional area 

L is the rod length 

X is the local coordinate system and longitudinal axis of the rod 

with origin at point 1, that is Xj^ ■ 0; X 2 • L. 

The equivalent mass at rod end 2 is calculated by f-king the first 
moment about rod end 1. 
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- I P ^ (A2-A^)) X dx 

Jo 

from which 

^ (A^+2 A^) 

and ■ M - M 2 

g— (2Aj^ + A2). 

The kinetic energy of the rod elemeut is expressed as 


[ 2 ] 

[31 


T 


1 





/x2<^>_ 


[41 


'ie kernel matrix of £q [4] is the mass matrix that represents the 
distributed mass of the rod. 
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Subroutine M1A2 calculates a consistent mass matrix for an axial 
rod element with unrestrained boundaries. The mass matrix Is In 
the local coordinate system of the rod. The elements of the mass 
matrix represent the distributed mass properties of the rod. 

These matrix elements are calculated by assuming the displacement 
between the rod ends to be a linear function of the displacement 
at the rod ends. The rod may be linearly tapered or uniform. 

DESCRIPTION OF TECHNIQUE 

The replacement of the distributed mass of a rod by a mass matrix 
is obtained using a klne''ic energy approach as follows. 

For small deflections (5) along the longitudinal axis (x) of the 
rod shown In the sketch, the kinetic energy Is defined by 


T 


1 

2 



pA(x) 6^ (x.t) dx 


U] 


where 

p is the mass density 
A Is the cross-sectional area 

<5 Is the time rate of change of displacement along the rod x-axls, 
referred to aa longitudinal velocity In the paper 


t is time 


X Is the local coordinate system and longitudinal axle of the r> d 
with or:'3ln at point 1, that is, - 0; X2 - L (rod length). 



Rod Element 
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To integrate Eq [1] a linear displacement function will be 
assumed between points 1 and 2 in terms of the displacements of 
points 1 and 2, that is. 


Hx) 


z 


(*2- ‘l>- 


Similarly, the longitudinal velocity is given by 


«(x,t) 


6l(t) + (6^(0 - 6^(t)) 


j- [(L-x) xj 


«^(t) 




[21 


The cross**sectionaI area of the rod is assumed to vary linearly, that 
is, 

V(x) » i (Aj-A^). [31 


Substituting Eq [2] and [3] into [1] gives the kinetic energy as 



r* 1 


“ “ 


— 

\(t) 

1 • • -T 

3Ai + A^ Aj^ + Aj 


6j(t) 

L(c) 

* I § 

Ai + A^ A^ + 3A2 


62 <t) 


The kernel matrix of £q [4] is the mass matrix. 


[4] 


The total mass properties of the rod may be calculated from the 
following triple matrix product. 


M P°' 


'l l‘ 


3A, ^ A_ A, + A« 


’i o' 


pL 



12 12 




12 






P° 1° 

_ m 


0 L 


Aj^ + A^ ^1 ^^2 


1 L 


w>^ere 

M is the mass of the rod 

is the first moment about Xj^ = 0 
is the moment of inertia about = 0. 


' 3 ] 
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Expanding Che triple matrix product gives 


M - p (Aj^+Aj) L/2 

[5a] 

P° - p (A^+Aj) L^/6 

I5bJ 

1° - p (A^+ 3A2> L^/12 

I5c] 

The center of gravity is calculated from 

X - P°/M 

eg 

(A^+ 2 A 2 ) L 
(Aj^+Aj) 3 

[5d] 
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Subroutine MlBl calculates a lumped mass matrix for a bending 
beam element with unrestrained boundaries. The mass matrix is in 
the local coordinate system of the beam. The elements of the mass 
matrix represent the distributed mass properties of the beam. 

These matrix elements are calculated by lumping the beam's total 
mass to the beam end points using static equivalent forces. The 
beam may be linearly capered or uniform. 

DESCRIPTION OF TECHNIQUE 

The replacement of the distributed mass of a beam by a mass matrix 
is obtained by lumping the beam's total mass Co the beam end points. 
This lumping is equivalent to static beaming. The beam is illus- 
trated in the sketch in which the rod cross-sectional area varies 
linearly, that is, A(x) “ A + x/L (A2-A ) . 



The tocal mass of the beam Is 

M = p L(A^+A2>/2 [1] 

where 

p is the mass density 

A is the cross-sectional area 

L is the beam length 

X is the local coordinate system and longitudinal axis of the 

beam with origin at point 1, that is, x^^ = 0; ^^2 ” 
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Th« equivalent mass at bevi end 2 Is calculated by taking the 
first moment about beam end 1. 

^2^ = I p Ax dx 
^o 

» I P (^1 L 

'o 

~6~ ^^2^ 

- pL (A^+ 2A2)/6 
From 

» M 

= pL (2 Aj^+A2)/6 

An attempt at calculating the inertias was made by taking the 
second moment about beam end 1. 


[21 


[3] 


II + I 2 ^2 “ 1 *** 

- ML^/3 

II + I 2 • (M^-ZM^) L^/3 
For a uniform beam, « M 2 ■ M/2 
+ I 2 » - ML^/6 


which is impossible. This says that lumping can never give the cor- 
rect inertia values. Therefore, arbitrarily assume 


= iA^ L^/24 


[A] 
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and 

- pAj L^/24 


The kinetic energy of the beam element is e3q>ressed as 


T 


y [same as coluam] 



The kernel matrix of [6] is the mass matrix that represents the 
distributed mass of the beam. 


[51 


16] 
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SubrouClne M1B2 calculates a consistent aass matrix for a bending 
beam element with unrestrained boundaries. The mass matrix Is In 
the local coordinate systoi of the beam. The elements of the mass 
matrix represent the distributed mass properties of the beam. 

The matrix elements are calculated by assuming the dlsplacanent 
between the rod ends to be a cubic function uf the displacement at 
the beam ends. The beam may be linearly tapered or uniform. 

DESCRIPTION OF TECHNIQUE 

The replacement of the distributed mass of a beam by a mass 
Biatrlx Is obtained using a kinetic energy approach as follows. 

For small deflections (6) normal to the longitudinal axis (x) of 
the beam shown in the sketch, the kinetic energy is defined by 

S *2 

p A(x) 6^ (x,t) dx [1] 

*1 

where 

p is the mass density 
A is the cross-sectional area 

6 is the time rate of change of displacement normal to the body 
x-axls, referred to as lateral velocity in this paper 

t is time 

x Is the local coordinate system and longitudinal axis of the 
beam with origin at point 1, that is, x. >0; x. ■> L (beam 
length) . 
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To integrate Eq [1] a cubic displacement function is assumed 
between points 1 and 2 in terms of the displacements of points 1 
and 2, that is. 


6(x) 


Ax^ + Bx^ + Cx + D 


x^ X 1] 


12 ] 


The angular displacement is obtained as the geometric derivative 
of the lateral displacement, that is 


e(x) 


_ a 5(x) 

3x 

t-3x^ -2x -10] 


A 


B 

C 

D 


[3] 


The coefficients A, B, C, D are determl4.«^ from Eq [2] and [3], 
using the displacements at the beam ends. 


'l" 


0 0 0 1 


A 



3 2 

L L 1 


B 


m 




CD 


0 0-10 


c 

CS 

<x> 


-3L^ -2L -1 0 


D 

1 




mm 


From which 


A 


r* "1 

B 


^2 

C 

= 1^1 


D 


mmi 


I4a] 
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where 




2/V 

l-3/I.^ 

0 

1 


-2/l" 

3/l' 

0 

0 


-1/l" 

2/L 

-1 

0 


-1/L^ 

1/T. 

0 

0 


[4b] 


Using Eq [4a] In [2] and taking the tine derivative gives the 
lateral velocity as 


fi(x,t) • [x^ x^ X 1] [♦! 


«^(t) 

ej^(t) 

e^Ct) 


The cross-sectional area of the beam is assumed to vary linearly, 
that is. 


[5] 


A(x) “ f (Aj-Aj^). 


[6] 


Substituting Eq [5] and [6] into [1] gives the kinetic energy as 




6 5 4 ,31 

XXX X ' 
X X X X 

4 J : 

X h X X 

X X X 1 


^ L 


M 


*i<') 


{si*r«e .» colannl 


T . 1. 
640 


— — 


T — 

2«*0A,+72A2 54A^+54A2 -OOA^+iAA^) 1 (14 Aj+12a^)L 


•j(t) 

,2Aj^ 240A2 -<l2Aj+14A2)i <l4Aj^+30A2>t 


•jU) 

U.\’n) (5 Aj-^3A2)L^ -(3Aj4-3A2>1^ 


yt) 

(3A^+5A2 C 







[7] 


The kernel matrix of Eq [7] is the mass matrix 
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Subroutine HlCl calculates a lumped aaas natrlx for a torsion rod 
elosent with unrestrained boundaries. The mass ipatrix is in the 
local coordinate systea of the rod. The elMmts of the oass 
matrix represent the distributed inertia properties of the rod. 

These matrix elei^tits are calculated by lumping the rod's total 
inertia to the rod end points using static equivalent forces. 

The rod may be linearly tapered or uniform. 

DESCRIPTION OF TECHNIQUE 

The replacoaent of the distributed inertia of a rod by a mass 
matrix is obtained by lumping the rods total inertia to the rod 
end points. This lumping is equivalent to static beaming. The 
rod is illustrated in the sketch in which the rod cross-sectlonal 
polar area moment of inertia varies linearly, that is, 

P(x) - + x/L 



The total inertia of the rod is 

I - pL (Pj+P2>/2 hi 


where 

p '*9 the mass density 

? the cross-sectional polar area moment of Inertia 
L Is the rod length 

X is the local coordinate system and longitudinal axis of the rod 
with origin at point 1, that is, *0; X 2 * L. 

The equivalent Inertia at rod end 2 is calculated by taking the 
f^rst moment about rod end 1. 
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'o 

»L 


pPx dx 


p (Pj^ + f (Pj-Pj)) * «ix 


Y ^ 

■S-k- rp •»• 2P ) 
6 ' 1 


from which 

T a fp +p ) 

2 6 1 2^ 


and 


I 


1 


I - 


I 


2 


- r«V‘'2>- 


The kinetic energy of 


T 


! ^ 2(01 


the rod element is expressed as 


- 


m m 



«xl<^> 

° ^2 


«x2<^> 

» m 




The kernel matrix of Eq [4] is the mass matrix that represents 
the distributed inertia of the rod. 


121 


[31 


[41 



Subroutlntt MlCl. calculates a consistent mass aatrlx for a torsion 
rod element with unrestrained boundaries. The mass matrix is in 
the local coordinate system at the rod. The elements of the mass 
matrix represent the distributed inertia properties of the rod. 
These matrix el«nents are calculated by assuming the displacemoat 
between the rod ends to be a linear function of the displacestent 
at the rod «ads. The rod may be linearly tapered or uniform. 

DESCRIPTION OF TECHNIQPE 

The replacement of the distributed Inertia of a rod by a mass 
obtained using- a kinetic energy approach as follows. 

rotations <6) along the longitudinal axis (x) of the 
in the sketch, the kinetic energy is defin^ by 


p P(x) (x,t) dx 


mass density 

cross-sectional polar area moment of inertia 

time rate of change of rotation along the rod x-axls, 
referred to as rotational velocity in this paper 

t Is time 

X is the local coordinate system and longitudinal axis of the 

rod with origin at point 1, that is, « 0; x^ * 1 (rod length) 


matrix is 

For small 
rod shown 



where 


p is the 
P is the 
6 is the 



Rod Elmnent 

To integrate Eq [1] a linear displacement function will be assumed 
between points 1 and 2 in terms of the rotations of points 1 and 2, 
that is, 

'j(x) “ L (®2"®1^' 

Similarly the rotational velocity is given by 
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0 (x,t) - e^(t) + I (ej(t) - e^(t) - ^ C(l-x) xJ - 




Oj(t)| 


The cross-sectional polar area moment of inertia of theorem is 
assisned to vary linearly, that is, 

P(x) -Pi+f<VV* 

Suhstituti g Eq [2] and [3] into [1] gives the kinetic energy as 


i-icvt) Vt))^ 

la 


ilU 


) I(L-X) x] dx 


0 j(t) 

8,(t) 




PL Px «2 ;,<t) 


(sym) KCt)! 


The kernel matrix of Eq [4] is the mass matrix. 

The total inertia of the rod can he calculated by assuming a rigid 
body mode of 


1 *2 

Substitution of [5] ipto [4] gives T * 1 0 » 


where 


3P,+P2 P,+P2 

I - [1 1] p^+P^ P1+3P2 


- 0 (p2+P2>/L2. 



H2A1 - 1/2 

Subroutin* M2A1 calculates a luaped aass natrlx for a nenbrane 
triangle plate eleiaent with unrestrained boundaries. The ness 
natrlx Is In the local coordlnata systen of the triangle plate. 

The eltf&ents of the nasa matrix represent the distributed mass 
properties of the triangle. 



The replacement of the distributed mass of a membrane triangle 
plate by a mass natrlx Is obtained by lumping the triangle's 
total mass to the triangle vertices. The triangle is Illustrated 
In the sketch with the degrees of freedom shown. 



The total mass of the triangle is 

M - pt yj/2 111 

where 

0 Is the mass density 
t Is the plate thickness. 

The mass at each vertex for a translation degree of freedom is 
M/3. Because any inertia at a vertex will always be "heavy", 
arbitrarily assume this value co be M/3 also. 
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The kinetic energy of the membrane triangle plate element Is 
expressed as 



The kernel matrix of Eq [2] is the mass matrix. 



M2A2 


Subroutine M2A2 calculates a consistent aass natrix for a aeBbrane 
triangle plate eleaent with unrestrained boundaries. The nass 
■atrix is in the local coordinate systM of the triangle plate. 

The eleaents of the eass aatrix represent the distributed aass 
properties of the triangle. I^ese aatrix eleaoits are calculated 
by assuaing a quadratic displaceaent field. 

DESCRIPTIOM OF TECHMIOOE 

The replaceaent of the distributed aass of a anbrane triangle 
plate by a aass natrix is described in D.M. 169 Linear Strain 
Hertrane Triangle Element by V. A. Benfield and C. S. Bodley. 

The triangle is illustrated in the sketch with the degrees of 
freedoa shown. 


^3 



fi5 6 e ], [6 6 e (A 6 e ],1 . 

X y a 1 ‘ X y z 2 x y z 3J 
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Subroutine M2B1 calculates a luaped nass matrix for a banding 
triangle plate eleaent with unrestrained boundaries. The mass 
matrix is in the local coordinate system of the triangle plate. 

The elements of the mass matrix represent the distributed mass 
properties of the triangle. 

DESCRIPTION OF TECHHIQUE 

The replacement of the distributed mass of a bending triangle 
plate by a mass matrix is <H>tained by lumping the triangle's total 
mass to the triangle vertices. The triangle is illustrated in the 
sketch with the degrees of freedom shown. 




The total mass of the triangle is 

M - pt x^^!2 [1] 


where 

0 is the mass density, and 
t is the plate thickness. 

The mass at each vertex for a translation degree of freedom is M/3. 
Because any inert 'a at a vertex will always be "heavy”, arbitrarily 
assume this value to be M/3 also. 
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Th« kinetic wiergy of tlM bending triangle plate ele a ent is 
expressed as 



The kernel aatrix of Eq 12] is the aass aatrix. 



M2B2 


Subroutine M2B2 cnlculntos • consistent ness nstrix for e bending 
trisngle piste elennt with imrestrsined boonderies. The ness 
■strix is in the locsl eoordinste systen of the trisngle piste. 
The elenents of tlM ness nstrix represent the distributed ness 
properties of the trisngle. TlMse nstrix elenents ere cslculsted 
by sssuaing s cubic displscenent field. 

DESQtlPTIWI OF TECHHIQOE 

The replscenent of the distributed ness of s bmding trisngle 
piste by s nsss nstrix uses the technique (essentislly) described 
in Trianglulccr Elemente in Plate Bending by C. P. Bssely, Y. K. 
Oteung* B. M. Irons, end 0. C. Zienkiewics; AFFDL-TR'-66-80, 
Novesber 1966. The trisngle is illustreted in the sketch with 
the degrees of freedon shown. 




The order of the degrees of freedom is 

[ [6 0 e ]- 16 e 0 ]- 16 e e ],| . 
‘ z X yu * z X y'2 ‘ z x y'3j 
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Subroutins MASIA calculatas « finite el — »n t —ss — trlx for en 
exiel rod ele—n t «rlth onrestrelned bounderies. 

The use —trlx, else (x6, is in the globel coordinete directions. 

The globel coordinete order for eech ele— nt is (U« V» W) joint 1, 
then Joint 2 Hhere U. V, W ere trensletions. If the Euler engles 
ere zero et e Joint, U * 6^, V * 6^, V • 6^. The iess — trix 

— y be either lushed or consistent. 

This ness —trix is coeputed by first celculeting e — ss —trix in 
the locel coordi— te syst— for either e lusped — ss —trix or e 
consistent — ss —trix. Euler angles are th— used to transform 
the — ss —trlx fr— the local coordi— te syst— to the globel 
coordi— te directions. Direction cosi— s are not needed. 

DESCRlPTKMt OF TECHHIQUE 

The calculation of the mass —trix in the gl<d>al coordi— te direc- 
tio— is accomplish as follows. First a mass —trix is calculated 
in the local coordi— te syst— and can be giv— as [m^] “ f"!! 

1*21 * 2 ? 

for either the luag>ed mass —trix (reference subroutine MlAl) or 
the consistent — ss —trlx (reference sid>routine M1A2) . The off- 
diagonal ter— are zero for the lumped — ss —trix. 

A sketch of the rod is given for reference as 



z 
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The deflections In the local systea are related to the deflections 
in the global c^rdlnate directions by 



where 


[y] is a direction cosine matrix, size 3x3, 

[£.] is an Euler angle transformation matrix (reference subroutine 
EULER) at joirt i, size 3x3, and 

{hg} = [U^ U 2 W^). U, V, W are translations. 

Substituting Sq. [2] into [1] gives the kinetic energy using global 
coordinates as 


T 




Iv]" 

[y] 

Hi) 

: ”12 
1 

“ *“ * 

[Ej]^ 

ty]^ 

ly] 

[E2] 

(Ej]^ 

[yf 

[y] 

lE^l 

1 ">22 

[£2]^ 

[y]^ 

[y] 

tEj] 


(h,) 


PI 
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The Unetic energy using local coordinates is 


1 T 

T - j [sane as coluan] 


"ll I "l2 

*21 i -22 “1 


»1 


yi 


xl 


y2 

*.2 


where [1] is a unity natrlx, size 3x3. 


til 


The deflections In the local system are related to the deflections 
in the global coordinate directions by 


xl 


yi 


zl 


^x2 

^2 


z2 


• 

* 

[y 1 lEj^l 

0 

0 

[y 1 [Ejl 


where 


{h^} 


121 


[y] is a direction cosine matrix, size 3x3, 

[E,l is an Euler angle transformation matrix (reference subroutine 
EULER) at Joint i, size 3x3, and 

{h^) * [Uj^ Vj^ Uj Vj W 2 I. U, V, W are translations. 

Substituting Bq. [2] into [1] gives the kinetic energy using global 
coordinates as 


1 * T 

"11 

[yI'^ 

tYl 

lE^l 

! “l2 

[E^l^ 

lYl’^ 

IyI IE 2 I 

T - f {h 


• ~ » 

M « 

mm mm ^ m 

1 

m M mm m 


> w w *■ 



®21 -^2^^ 


[Y] 

(Ej) 

B 

to 

to 

[£ 2 !’’ 

[y]^ 

Iy] [E 2 ] 


I 


[31 




For luapod aass, kernel aetrlx in Eq [4] is the 

desired mss Mtrix in the globel coordinete directimis. If the 
Euler sngles ere sero at both Joints, [E^] ■ {!]. 
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Sud>routi&e >US1B calculatas a fiaite el«nant mass matrix for a 
comblnad axial-*tor8lon>bending bar alemant with unreatrainad 
boundaries. 

The aasa matrix* size 12x12 > is in the global coordinate direc- 
tions. The global coordinate order for each element is (U, V, 

W* P, Q* R) Joint 1, then Joint 2 where V* V, W are translations 
and P* Q* R are rotations. If the Euler angles are sero at a 
joint* then U*6 * V»6 * W«6 , P«»6 * Q»6 * R»9 . 

X Y* Z X Y Z 

This mass matrix is computed by first calculating a mass matrix 
in the local coordinate system for either a lumped mass matrix 
or a consistent mass matrix. A direction cosine matrix is then 
used to transform the mass matrix from the local coordinate sys- 
tem to the global coordinate directions. 

DESCRIPTION OF TECHNIQUE 

The calculation of the mass matrix in the global coordinate direc- 
tions is accomplished as follows. First a mass matrix is calculated 
in the local coordinate system for either limtped mass or consistent 
mass using uncoupled axial* torsion* and bending subroutines 
listed. 



lumped 

consistent 

axial 

MlAl 

HU2 

bending 

MlBl 

M1B2 

torsion 

MlCl 

M1C2 


A sketch of the bar is given for reference as 



The kinetic energy using local coordinates is 
I - I (1.^1 


m 
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where 

T ' II I 

L ^ xl x2 I xl x2 I yl y2 si z2 I xl z2 yl y2^ 


and 



refer to terms In the uncoupled axial, torsion, bending 
mass matrices. For lumped mass, the off-diagonal terms are zero. 

The above deflections in the local coordinate directions are related 
to deflections in the global coordinate directions by 

{hj^} - IyI {hg} [2] 

where (y 1 is a direction cosine matrix (reference subroutine DCOSIB) 
including Euler angles, size 12x12, and {h^}^ ■ 

^2^2^2^2^2^2 ^ ^ translations and P, Q, R are rotations. 
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Substituting Eq [2] into [1] givss ths kinstic snergy using 
global coordinates as 

T - i {hg}’' [mg] (hg) [3] 

T 

where [Wg] • [y] [i^] [y] ie the desired mass matrix in the 

global coordinate uirections. 

Even though the local lumped mass matrix has only diagonal terms, 
the triple matrix product using direction cosines is needed because 

'u '* ‘as “•* '22 ’* ‘m- 
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Subroutln* MAS2 calculates a finite elotcnt mass matrix for a 
combined mambrana-bandlng triangle plate alament with unrestrained 
boundaries. 

The mass matrix, slse 18x18, Is In the global coordinate directions. 
The global coordinate order for each element Is (U, V, W, P, Q, R) 
joint 1, then joints 2, 3 where U, V, W are translations and P, 

Q, R are rotations. If the Euler angles are zero at a joint, then 
U-6 , V-6 , W-« , P-8 , Q-8 , R-e . 

X T X X T Z 

This mass matrix Is computed by first calculating a mass matrix 
in the local coordinate system for either a lung>ed mass matrix 
or a consistent mass matrix. A direction cosine matrix Is then 
used to transform the mass matrix from the local coordinate sys- 
tem to the global coordinate directions. 

DESCRIPTION OF TECHNIQUE 

The calculation of the mass matrix in the global coordinate directions 
is accomplished as follows. First a mass matrix is calculated in 
the local coordinate system for either lumped mass or consi& *:ent 
mass using uncoupled membrane and bending subroutines listed. 



lumped 

consistent 

membrane 

banding 

M2A.1 

M2B1 

M2A2 

M2B2 


A sketch of the triangle plate is given for r^ ference as 
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The kinetic energy using local coordinates is 


where 


<\> - ‘,1 

®zl 




and 




’iv 

1 

1 

1 

Kl - 

mem 

1 

» i 

e 


0 


x2 ^y2 ®z2 

*x3 *y3 *.3j 

z2 ®x2 ®y2 

«.3 *x3 «y3> 


(«,] 

ben 


[11 


The above deflections in Che local coordinate directions are 
related to deflections in the global coordinate directions by 

{hj^} - lYl {h^} [21 

where [y] is e direction cosine matrix (reference subroutine 
DC<^S2) including Euler angles, size 18x18, and {h^}l ■ 

^^1 '^l ”l ^1 ^1 ^1 ^2 ^^2 ”2 ^2 ^^2 *^2^' translations 

and P, Q, R are rotations. 

SubstlCutlng Eq [2] into [1] gives the kinetic energy using 
global coordinates as 

T - I (m^l {h^} [31 

T 

where [m^] ~ [y] [m, ] [y] is Che desired mass matrix in the 

, 7 . 1 obal coordinate directions. 

Because the local lumped mass matrix is diagonal and because 
translation terms are equal for membrane and bending, and rotation 
terms are equal for membrane and bending, the triple matrix product 
using direction cosines is not needed. A simple reordering of 
diagonal terms is sufficient. 



M\S3 

Subroutine MAS2 calculates a finite element mass matrix for a 
combined membrane-bending quadrilateral plate element with un- 
restrained boundaries. 

The mass matrix, size 24x24, Is In the global coordinate directions. 

The global coordinate order for each element Is (U, V, W, P, Q, K; 
joint 1, then joints 2, 3, 4 where U, V, W, are translations and 
P, Q, R are rotations. If the Euler angles are zero at a joint, 

then U*6 , V-6 , W-<S , P-6 , Q-6 , R-0 . 

X Y Z X Y Z 

This mass matrix Is calculated by taking the average overlap of 
four triangles show.. In the sketch. Either a lumped mass matrix 
or a consistent mass matrix Is calculated. Subroutin'. MAS2 Is 
used for the calculation of the mass matrix for the triangular 
places. 


3 



Z 



(^lAD 


Subroutine QUAD calculates (on option) finite el«ent: (1) nass 

■atrices; (2) stiffness eatrices (saaw as global load transforaation 
■atrices), and (3) vectc. s to locate the DOF (degrees of freedom) 
of the aatrices in the global DOF. for coad>ined aesbrane-bending 
quadrilateral plate elements. The above aatrices and vectors are 
written on disk units and constitute the output from this sub- 
routine. All aatrices are in dense progr awning logic. 

Each mass and stiffness aatrix, size 24x24. is in the global 
coordinate directions. The global coordinate order for each ele- 
ment is (U, V. U, P. Q, R) Joint 1, then Joints 2. 3. and 4 
where U, V, U are translations and P, Q. R are rotations. If 

the Euler angles are zero at a Joint, then U>5 . V>6 . W»6 . P^d . 

X Y Z X 

Q»9 . R*0 . 

Y Z 

Each global load ti : . raatlon catrix. size 24x24. relates loads 

at quadrilateral . ..ce-i in the global coordinate directions to 
deflections in the global coordinate directions. The row order 
in this matrix is (P^. P^. P^. H.. M^. »^) Joint i. then Joints 2, 3 

and 4 vheie P is force end M Is raomente 

Each location vector (IVEC) locates the DOF of each finite element 
in the global DOF. For exai^>le« 1VEC(6)«834 places element DOF 6 
into global DOF 834. 1VEC(3)^0 OaJ.ts element DOF 3 from global 

JOF. This constrains element DOF 3 to zero motion. 

The above matrices are calculated by using Joint data and element 
data. The Joint data is obtaxned from three matrices input to 
this subroutine: (1) Joint global Xp Y» Z locations; (2) Joint 

global DOF numbers; and (3) Joint Euler angles. 

The element data» read in this subrv «tinep is: (1) options for 

tnassp stiffness; (2) element material properties; and (3) element 
joint numbers. 


Each mass matArix is calculated by transfer to subroutine MAS3. 

Exch stiffness matrix is calculated by transfer to subroutine STF3. 
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Subroutine STFIA calculates a finite element; (1) stiffness 
matrix (same as global load transformation matrix) ; (2) local 
load transformation matrix; and (3) on option, stress transforma- 
tion matrix for an axial rod element with unrestrained boundaries. 

The stiffness matrix, sixe 6x6, is in the global coordinate 
directions. The global coordinate order for each element is 
(U,V,U) joint 1, then joint 2 where U,V,U are translations. If 

the Euler angles are zero at a joint, then , V«6 , U”6 . 

X T Z 

The global load transformation matrix, size 6x6, relates loads 
at the rod ends in the global coordinate directions to deflections 
in the global coordinate directions. The row order in this 
matrix is (Py, P^, P^') joint 1, then joint 2 where P is force. 

The local load transformation matrix, size 2x6, relates loads 
at the rod ends in the local coordinate system to deflections in 
the global coordinate directions. The row order in this matrix 
is P^j^, P ^2 where P^ is axial force. 

The stress transformation matrix, size 2x6, relates stresses at 
the rod ends in the local coordinate system to deflections in the 
global coordinate directions. The row order in this matrix is 
^xl* ^x2 ^ normal stress. 

These matrices are computed by first calculating a stiffness 
matrix and stress transformation matrix in the local coordinate 
system. A direction cosine nmtrix is then used to transform 
the stiffness matrix and, on option, the stress transformation 
matrix from the local coordinate system to the global coordinate 
directions. 


DE SCRIPTION OF TECHNIQUE 


The calcilation of the stiffness matrix, load transformation matrix, 
and stress transformation in the global coordinate dlrectic is is 
accomplished as follows. First a stiffness matrix is calculated 
in the local coordinate system and can be given as [k^] = 

(reference subroutine KlAl) . 


Ic k. 

11 12 


Ic Ic 
21 22 


A sketch of the rod is given for reference as 



X 
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The strain energy using local coordinates is 




The deflections in Che local systen are related to the deflections 
in Che global coordinate directions by 


'^xl 

- <V ■ 

x2 ^ 


[21 


where [y] Is a direction cosine matrix (reference subroutine 
DCOSIA) including Euler angles, size 2x6, and 

{h_}^ = [U- V W U_ V- W ]- U, V, W are translations. Substituting 
G i. X X 2 Z 

Eq [2] into [1] gives 


u . i (hg)- ;k„l (hg) 13] 

T 

where [k^] = [7] [k^] [y] is the stiffness matrix in global 

coordinate directions. 

The loads in the global coordinate directions can be calculated from 
Eq [3] as 
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Thus [k-] is also a global load transforaation awtrix giving loads 

in the global coordinate directions to deflections in the global 
coordinate directions. 

The loads in the local coordinate directions can be calculated 
from Eq [1] as 


“ Ffhjr • 

Substituting Eq [2] gives 

- tTLj {hg} [6] 

where [TL] « [li^] [y 1 is the local load transformation matrix 

giving the loads in local coordinate directions to deflections in 
the global coordinate directions. 

A stress transformation matrix relating stresses in the local 
coordinate directi ns to deflections in the local coordinate 
directions is first calculated (reference subroutine KlAl) , that 

is, 

ISl> = (TSj^] {hj^} . [7] 

On option, the stress transformation matrix relating stresses in 
the local coordinate directions to deflections in the global 
coordinate directions is calculated. Substituting Eq [2] into 
[71 gives 

is^) = [IS] {hg} 
where 

(TSl - (TSj^l (y1 . 
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Subroutine STFIB calculates a finite element: (1) stiffness 
matrix (same as global load transformation matrix) ; (2) local 
load transformation matrix; and (3) on optiont stress transforutlon 
matrix for a combined axial-torsion-bendlng bar element with un- 
restrained boundaries. 

The stiffness matrix, size 12x12. is in the global coordinate 
directions. The global coordinate order for each element is (U, V. 
W, P, Q. R) joint 1; then Joint 2 «rhere U, V. U are translations 
and P» Q, R are rotations. If the Euler angles are zero at a 

joint, then U»6 , V-6 , W*6 , P=6 , Q=0 , R=8 . 

X T Z X Y’ Z 

The global load transformation matrix, size 12x12, relates loads 
at the bar ends in the global coordinate directions to deflections 
in the global coordinate directions. The row order in this matrix 
is (Py, Py Py, Mp, Mq, joint 1; then joint 2 wliere P is force 

and M is moment. 


The local load transformation matrix, size 12x12, relates loads 
at the bar ends in the local coordinate system to deflections in 
the global coordinate directions. The row order in this matrix 


is P 


M 


xl’ ‘x2* “xl* 

P is force and M is moment 


\2- V* 


n 


y2’ “zl* ”z2» "zl’ *^zl* "yl* "y2* 


M 


M 


where 


The stress transformation matrix, size 12x12, relates stresses 
at the bar ends In the local coordi ate system to deflections in 
the global coordinate directions. The row order in this matrix is 


xl 


x2 


M * r M * r 
xl 1 x2 2 


P P 

!zl 1X2 

A. ’ A„ 


\l * %1 

^zl 


^z2 * %2 
^z2 


zl 

where 


z2 


M * c , 
yl zl 


\2 * ^z2 

^2 


P is force, 

M is moment, 

A is cross-sectional area 

r is distance from torsion axis to outer fiber 

ic c.io**.3-sectiou Saint Venant's torsion constant JG 
c is distance fi^om bending m^utral plane to outer fiber an-* 
1 is area moment of inertia about local axis. 
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These eatrlces are c<»puted by first calculating a stiffness 
matrix and stress transformation isatrix in tne local coordinate 
system. A direction cosine matrix is then used to transform the 
stiffness matrix and, on option, the stress transformation matrix 
from the local coordinate system to the global coordinate directions. 

DESCRIPTION OF TECHNIQPE 

The calculation of the stiffness matrix, load transformation matrix, 
and stress transformation in the global coordinate directions Is 
accoq>lished as follows. First a stiffness matrix is calculated in 
the local coordinate system using uncoupled axial, torsion, and 
bending subroutines listed. 



Subroutine 

axial 

KlAl 

bending 

KlBl 

torsion 

Kiel 


A sketch of the bar is given for reference as 



Strain energy using local coordinates is 

U . i ,h^) (k^l (hj^> 

where 


{h, } e , 9 . : -i , 6 . e , 9 . : ^ ^ ^ 1 ® 

L xl x2 . xl x2 ; yl y2 zl z2 ; zl z2 yl y2 


[ 1 ] 
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and 



stiffness matrices. 

Deflections in the local system are related to deflections in the 
global coordinate directions by 


{\} - [y 1 {h^} C21 

where [y3 is a di rection cosine matrix (reference subroutine DCOSIB) 
including Euler angles, size 12x12, and 

{h^}^ « [U^ U 2 '^2 ’^2 ^2 ^2^* U, V, W are transla- 

tions and P, Q, R are rotations. 

Substituting Eq [2] into [1] gives 

U . i (h,)" [k^l (h,) t31 

where [k„] * [v]^ (k. ] [y] is the stiffness matrix in global 
G L 

coordinate directions. 
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Loads in global coordinate directions can be calculated from 
Eq [3] as 

■ Ith^r ‘“g' ‘"c* “> 

Thus, [kg] Is also a global load transformation matrix giving 

loads in the global coordinate directions to deflections in the 
global coordinate directions. 

Loads in local coordinate directions can be calculated from Eq [1] 
as 

<»L> ■ TTiljT <N.> 

Substituting Eq [2] gives 

{pj^} - [TL] {hg} [6] 

Where [TLl = [kj^] [vl is the local load transfonoation matrix giving 

the loads in local coordinate directions to deflections In the 
global coordinate directions. 

A stress transformation matrix relating stresses In local coordinate 
directions to deflections In local coordinate directions Is first 
calculated (reference subroutines KlAl, KlBl, KIC) , that Is, 

{SlI * [TSj^] {h^} . [7] 

On option, the stress transformation matrix relating stresses in 
local coordinate directions to deflections In global coordinate 
directions is calculated. Substituting Eq [2] into [7] gives 

{sj^} = [TS]fhg} 
where 

ITS) - [TS^] Iy). 
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Subroutine STF2 calculates a finite element: (1) stiffness 

matrix (same as global load transformation matrix) ; (2) local lead 
transformation matrix; and (3) on option, stress transformation 
matrix for a combined membrane-bending triangle plate element 
with unrestrained boundaries. 

The stiffness matrix, size 18x18, is in the global coordinate 
directions. The global coordinate order for each element is (U, 

V, W, P 0, R) joint 1; then joints 2 and 3 where U, V, W are 
translations and P, Q, R are rotations. If the Euler angles are 

zero at a joint, then U“6 , V*6 , W»6 , P*8 , Q*6 , R»6 . 

X Y Z X Y Z 

The global load transformation matrix, size 18x18, relates loads 
at the triangle vertices in global coordinate directions to de- 
flections in global coordinate directions. The row order in this 
matrix is (P^, P^, P^^, Mp, Mq, Mj^) joint 1; then joints 2 and 3 

where P is force and M is moment. 

The local load transformation matrix, size 18x18, relates loads 
at the triangle vertices in the local coordinate system to deflec- 
tions in global coordinate directions. The row order in this 
matrix is (P , P , M ) joint 1; then joints 2 and 3, next ^P , M , 

My) joint 1; then joints 2 and 3 where P is force and M is moment. 

The stress transformation matrix, size 18x18, relates stresses at 

the triangle vertices in the local coordinate system to deflections 

in global coordinate directions. The row order in this matrix is 

(j , c , T ) at 2 = + l/2t. for joint 1; then joints 2 and 3. 

X -V xy ben '' 

Then the same data at Z = - o is normal stress and x is 

shear stress. 

These matrices are computed by first calculating a stiffness matrix 
and stress transformation matrix in the local coordinate system. 

A direction cosint matrix is then used to transform the stiffness 
matrix and, on opt. on, the stress transformation matrix from the 
local coordinate system to the global coordinate directions. 

DESCRIPTION OF TECHNIQUE 

The calculation of the stiffness matrix, load transformation 
matrix, and stress transformation in global coordinate directions 
is accomplished as follows. First a stiffness matrix is calculated 
in the local coordinate system using the uncoupled membrane and 
bending subroutines listed. 



Subroucine 

membrane 

bending 

K2A1 

K2B1 



A sketch of the triangle platw Is given for reference as 
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Z 

The strain energy ’ising local coordinates is 


U - I [kj^] {1^} [1] 


where 


*yl ®zl ^x2 ^y2 ®z2 ^x2 ^y3 ®z3 


6 , 0 , 6 . 5 - 0 - 0 . 6 , 6,6 1 

zl xl yl z2 x2 y2 z3 x3 y 


and 



I 

I 

I 

I 

r 

I 

1 

I 




STF2 - 3/4 


Deflections in the local systes are related to the deflections 
in global coordinate directions by 

- {hj^} - lYl {hg} [2] 

where [y] Is a direction cosine matrix (re. . ''ence subroutine DC0S2) 
Including Euler angles, size 18x18, and 

(hc)^ - 1«1 *1 «1 0 i h "a *2 "a ‘'2 ‘*2 *2>- 

U, V, W are translations; P, Q, R are rotations. 

Substituting Eq [2] into [1] gives 

U » I Ik^J {hg} [31 

T 

where [k^] ~ [y] [k^l [y1 is Che stiffness matrix in global 
coordinate dire^'Clons. 

Loads in the global coordinate directions can be calculated from 
Eq [3] as 

“ ~{'h3 “ 

Thus, [k„] is so a global load transformation isatrix giving loads 

in global coordinate directions to deflections in global coordinate 
directions. 

Loads in the local coordinate directions can be calculated from 
Eq [1] as 

-IV' 

Substituting Eq [2] gives 

{Pj^} - [TL] {hg} [6] 

where [TL] » [k^^] [y] is the local load transformation matrix giving 

the loads in local coordxnate directions to deflections in global 
coordinate directions. 
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A stress transfomstion matrix relating stresses in the Iccal 
coordinate directions to deflections in the local coo linate 
directions is first calculated (reference subroutines K2A1, K2B1) , 
that is, 

{Sj^} - (TSj^l {hj^} . [7] 

On option, the stress transformation matrix relating stresses in 
local coordinate directions to deflections in global coordinate 
directions if. calculated. Substituting Eq ( 2] into [7] gives 

{Sj^} - (TSHhg} 

where 

(TS] - [y1. 



ST?3 

Subroutine STF3 calculates a finite element stiffness matrix 
(same as global load transformation matrix), for a coniblned 
membrane-bending quadrilateral plate element with unrestrained 
boundaries. 

The stiffness matrix, size 14x24, is in global coordinate direc- 
tions. The global coordinate order for each element is (U, V, W, 

P, Q, R) joint 1; then joints 2, 3 and 4 where U, V, W are transla- 
tions and P, Q, R are rotations. If the Euler angles are zero at 
a joint, then U“6 , V«: , W»6 , P»0 , 0*0 , R“0 . 

X y’ Z X Y z 

Each global load transformation matrix, size 24x24, relates loads 
at the quadrilateral vertices in global coordinate directions to 
deflections in global coordinate directions. The row order It. 
this matrix is (P^, P^, P^, Mp, M^, Mj^) joint 1; then joints 2, 

3 and 4 where P is force and M is moment. 

This stiffness matrix is calculated by taking the average overlap 
of four triangles, shown in the sketch. Subroutine STF2 is used 
for calculation the. stiffness nutrlx for the triangular pastes. 


3 





2 
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Subroutine TRMGL calculates (on option) finite element: (1) mab 
matrices, (2) stiffness matrices ^ame as global load transformation 
matrices); (3) local load transformation matrices; (4) stress trans- 
formation matrices; and (5) vectors to locate the DOF (degrees of 
freedom) of these matrices In the global DOF, for combined mend>rane- 
bendlng triangle plate elements. The<^e matrices and vectors are 
written on disk units and constitute t .e output from this 
routine. All matrices are In dense progranning logic. 

Each mass and stiffness matrix, size 18x18, Is In global cocrdinaro 
directions, xhe global coordinate order for each elcmer.t is (L, 

V, W, P, Q, F) joint 1; then joints 2 and 3 where U, V, gre 
translations and P, Q, R are rotations. If the Euler i>... 'as arc 

zero at a joint, then U=6 , V«6 , W«(S , P«»3 , Q=6 , R=6 . 

X Y Z X Y Z 

Each global load transformation matrix, size 18x18, relates loads 
at the triangle vertices In global coordli ate dlreccions to 
deflections In global coordinate directions. The row order in 
this macrlx Is (Py, P^, P^, M^, M^, M^^) joint 1; then joints 2 

and 3 where P Is force and M Is moment. 

Each local lo..- ’ transformation matrix, size 18xlo, relates loans 
at the triangle vertices In the local coordinate system to deflec- 
tions in global coordinate directions. The row order In this 

matrix is (P , P , M ) joint 1; then joints 2 and 3 ne?. (P , M , 
x y z z X 



Each stress transformation matrix, size 18x18, relates stresses 
at the triangle vertices in the local coordinate system to deflec- 
tions ■'n glo’jal coordinate directions. The row order in this 

matrix is (o , < , x ; at 7 + l/2t. for joint 1; then joints 

:i y xy ben ■’ ■’ 

2 and 3; then the same data at Z * -l/2t, . o is normal stress 

ben 

and T is shear stress. 

Each location sector (IVIC) locates the DOF of each finite ele- 
ment in the global DOF. For example, IVZC(6)“834 places element 
DOF 6 into global DOF 834. IVEC(3)=*0 . cs element DOF 3 from 

global DOF. This constrains element DOF 3 to zero motion. 

The ctaove matrices are calculated by uaing joint data and eleir.cn 
data. The joint data, obtained from three matrices input to .nia 
subroutine, are il) joint global X, Y, Z i. cations, (2) joint 
global DOF numbers, and (3) joint Euler angles. 

The element data read in this sc routine lb (1) options for mas". , 
stiffness, local load transformations, stress cransformatlori~ ; 

(2) element material properties; and (3) element joint numbers. 



TBHGL - Vt 

Bach mam» ipttcix !• ,c*lc tt l« f d by tr«wf«r to oubroutliie MAS2. 

Boe2 - .If f aoOT •aurlx, loote and stroM trMsforMtion mtrlx is 
csica.4ttad by trsuf«r to sdbxo^ino STF2. 



